function [M,param] = dasmtx(SIG,x,z,varargin)

%DASMTX   Delay-and-sum matrix
%   M = DASMTX(SIG,X,Z,DELAYS,PARAM) returns the numel(X)-by-numel(SIG)
%   delay-and-sum DAS matrix. The matrix M can be used to beamform SIG (RF
%   or I/Q signals) at the points specified by X and Z.
%
%   Try it: enter "dasmtx" in the command window for an example.
%
%   Because the signals in SIG are NOT required (only its size is needed)
%   to create M, the following syntax is recommended:
%       M = DASMTX(size(SIG),X,Z,DELAYS,PARAM)
%   !IMPORTANT! -- With this syntax, use M = DASMTX(1i*size(SIG),...) to
%   return a complex DAS matrix for I/Q data.
%
%   NOTE: SIG must be a 2-D or 3-D array (the 3rd dimension is ignored).
%         Or, size(SIG) must be a two- or three-component vector (the 3rd
%         component is ignored). The first dimension (i.e. each column)
%         corresponds to a single RF or I/Q signal over (fast-) time, with
%         COLUMN #k corresponding to ELEMENT #k.
%
%   DASMTX returns the same results as DAS.
%       1) using DAS:
%           bfSIG = das(SIG,x,z,delays,param,method);
%       2) using DASMTX:
%           M = dasmtx(size(SIG),x,z,delays,param,method);
%           bfSIG = M*SIG(:);
%           bfSIG = reshape(bfSIG,size(x));
%
%   DELAYS are the transmit time delays (in s). The number of elements in
%   DELAYS must be the number of elements in the array (which is equal to
%   size(SIG,2)). If a sub-aperture was used during transmission, use
%   DELAYS(i) = NaN if element #i of the linear array was off.
%
%   PARAM is a structure that contains the parameter values required for
%   DAS beamforming (see below for details).
%
%   M is a large sparse matrix. Computing M can be much more cost-effective
%   than using DAS if you need to beamform several SIG matrices, because
%   M needs to be determined only once.
%
%   Let us consider that a series SIG{1}, SIG{2} ... SIG{N} of ultrasound
%   matrices have been generated by sending similar wavefronts with the
%   same array. These signals SIG{i} are stacked in a 3D array sig3D so
%   that sig3D(:,:,i) = SIG{i}. To beamform these data with a delay-and-sum
%   approach, the following can be used:
%       M = dastmtx([size(sig3D,1) size(sig3D,2)],x,z,delays,param);
%       % (or M = dastmtx(1i*[size(sig3D,1) size(sig3D,2)],...) for I/Q data)
%       bfSIG3D = M*reshape(sig3D,[],size(sig3D,3));
%       bfSIG3D = reshape(bfSIG3D,size(x,1),size(x,2),[]);
%
%   You can also consider saving the DAS matrix M in a MAT file and loading
%   it when needed. The previous syntax is much faster than:
%       for k = 1:N
%           bfSIG{k} = das(SIG{k},x,z,delays,param);
%       end
%
%   DASMTX(SIG,X,Z,PARAM) uses DELAYS = param.TXdelay.
%
%   DASMTX(...,METHOD) specifies the interpolation method. The available
%   methods are decribed in NOTE #3 below.
%
%   [M,PARAM] = DASMTX(...) also returns the structure PARAM with the
%   default values.
%
%   ---
%   NOTE #1: X- and Z-axes
%   The migrated signals are calculated at the points specified by (X,Z).
%   Conventional axes are used:
%   i) For a LINEAR array, the X-axis is PARALLEL to the transducer and
%      points from the first (leftmost) element to the last (rightmost)
%      element (X = 0 at the CENTER of the transducer). The Z-axis is
%      PERPENDICULAR to the transducer and points downward (Z = 0 at the
%      level of the transducer, Z increases as depth increases).
%   ii) For a CONVEX array, the X-axis is parallel to the chord and Z = 0
%       at the level of the chord.
%   ---
%   NOTE #2: DASMTX uses a standard delay-and-sum.
%   ---
%   NOTE #3: Interpolation methods
%   By default DASMTX uses a linear interpolation to generate the DAS
%   matrix. To specify the interpolation method, use DASMTX(...,METHOD),
%   with METHOD being:
%      'nearest'   - nearest neighbor interpolation
%      'linear'    - (default) linear interpolation
%      'quadratic' - quadratic interpolation
%      'lanczos3'  - 3-lobe Lanczos (windowed sinc) interpolation
%      '5points'   - 5-point least-squares parabolic interpolation
%      'lanczos5'  - 5-lobe Lanczos (windowed sinc) interpolation
%
%   The linear interpolation (it is a 2-point method) returns a matrix
%   twice denser than the nearest-neighbor interpolation. It is 3, 4, 5, 6
%   times denser for 'quadratic', 'lanczos3', '5points', 'lanczos5',
%   respectively (they are 3-to-6-point methods).
%   ---
%
%   PARAM is a structure that contains the following fields:
%   -------------------------------------------------------
%   1)  PARAM.fs: sampling frequency (in Hz, REQUIRED)
%   2)  PARAM.pitch: pitch of the transducer (in m, REQUIRED)
%   3)  PARAM.fc: center frequency (in Hz, REQUIRED for I/Q signals)
%   4)  PARAM.radius: radius of curvature (in m, default = Inf, linear array)
%   5)  PARAM.TXdelay: transmission delays (in s, required if DELAYS is not given)
%   6)  PARAM.c: longitudinal velocity (in m/s, default = 1540 m/s)
%   7)  PARAM.t0: start time for reception (in s, default = 0 s)
%
%   A note on the f-number
%   ----------------------
%   The f-number is defined by the ratio (depth)/(aperture size). A null
%   f-number (PARAM.fnumber = 0) means that the full aperture is used
%   during DAS-beamforming. This might be a suboptimal strategy since the
%   array elements have some directivity.
%
%   Use PARAM.fnumber = [] to obtain an "optimal" f-number, which is
%   estimated from the element directivity (and depends on fc, bandwidth,
%   element width):
%
%   8)  PARAM.fnumber: reception f-number (default = 0, i.e. full aperture)
%   9)  PARAM.width: element width (in m, REQUIRED if PARAM.fnumber = [])
%        or PARAM.kerf: kerf width (in m, REQUIRED if PARAM.fnumber = [])
%        note: width = pitch-kerf
%   10) PARAM.bandwidth: pulse-echo 6dB fractional bandwidth (in %)
%            The default is 75% (used only if PARAM.fnumber = []).
%
%   Advanced option for vector Doppler (Reception angle):
%   ---------------------------------------------------
%   11) PARAM.RXangle: reception angles (in rad, default = 0)
%       This option can be used for vector Doppler. Beamforming with at
%       least two (sufficiently different) reception angles enables
%       different Doppler directions and, in turn, vector Doppler.
%
%   Passive imaging
%   ---------------
%   12) PARAM.passive: must be true for passive imaging (i.e. no transmit).
%       The default is false.
%
%
%   REFERENCES:
%   ----------
%   1) If you use DAS or DASMTX, please cite:
%      V Perrot, M Polichetti, F Varray, D Garcia. So you think you can
%      DAS? A viewpoint on delay-and-sum beamforming. Ultrasonics 111,
%      106309. <a
%      href="matlab:web('https://www.biomecardio.com/publis/ultrasonics21.pdf')">PDF here</a>
%   2) If you use PARAM.RXangle for vector Doppler, please also cite:
%      Madiena C, Faurie J, Porée J, Garcia D. Color and vector flow
%      imaging in parallel ultrasound with sub-Nyquist sampling. IEEE Trans
%      Ultrason Ferroelectr Freq Control, 2018;65:795-802. <a
%      href="matlab:web('https://www.biomecardio.com/publis/ieeeuffc18a.pdf')">PDF here</a>
%
%
%   Example:
%   -------
%   %-- Generate RF signals using a phased-array transducer
%   % Phased-array @ 2.7 MHz:
%   param = getparam('P4-2v');
%   % TX time delays (80-degree-wide diverging wave)
%   dels = txdelay(param,0,80/180*pi);
%   % Scatterers' position:
%   xs = [(-1:0.5:1)*4e-2 zeros(1,5)];
%   zs = [ones(1,5)*6e-2 (2:2:10)*1e-2];
%   % Backscattering coefficient
%   BSC = [ones(1,9) 0];
%   % RF signals:
%   param.fs = 4*param.fc; % sampling frequency
%   RF = simus(xs,zs,BSC,dels,param);
%   % Plot the RF signals
%   subplot(121)
%   plot((RF(:,1:7:64)/max(RF(:))+(1:10)*2)',...
%      (0:size(RF,1)-1)/param.fs*1e6,'k')
%   set(gca,'XTick',(1:10)*2,'XTickLabel',int2str((1:7:64)'))
%   title('RF signals')
%   xlabel('Element number'), ylabel('time (\mus)')
%   xlim([0 22]), axis ij
%
%   %-- Demodulation and beamforming
%   % Demodulation
%   IQ = rf2iq(RF,param);
%   % Beamforming 256x128 80-degrees wide polar grid
%   [x,z] = impolgrid([256 128],9e-2,80/180*pi,param);
%   % DAS matrix and beamformed IQ
%   Mdas = dasmtx(1i*size(IQ),x,z,dels,param);
%   IQb = reshape(Mdas*IQ(:),size(x));
%   % Beamformed image
%   subplot(122)
%   pcolor(x*100,z*100,abs(IQb).^.5)
%   colormap(gray)
%   title('Gamma-compressed image')
%   xlabel('[cm]'), ylabel('[cm]')
%   shading interp, axis equal ij tight
%
%
%   This function is part of <a
%   href="matlab:web('https://www.biomecardio.com/MUST')">MUST</a> (Matlab UltraSound Toolbox).
%   MUST (c) 2020 Damien Garcia, LGPL-3.0-or-later
%
%   See also DAS, DASMTX3, SIMUS, TXDELAY, RF2IQ, GETPARAM.
%
%   -- Damien Garcia -- 2017/10, last update 2020/10/21
%   website: <a
%   href="matlab:web('http://www.biomecardio.com')">www.BiomeCardio.com</a>

%
%   PARAM can also contain the "elements" field:
%   -------------------------------------------
%   13) PARAM.elements: coordinates of the transducer elements (in m)
%       PARAM.elements must contain the x- (and optionally z-) coordinates
%       of the transducer elements. It must be a vector (or a matrix with 2
%       rows or 2 columns) corresponding to the x- (and optionally z-)
%       coordinates, respectively.
%       

if nargin==0
    if nargout>0
        [M,param] = RunTheExample;
    else
        RunTheExample;
    end
    return
end

%------------------------%
% CHECK THE INPUT SYNTAX %
%------------------------%

assert(nargin>3,'Not enough input arguments.')
assert(nargin<7,'Too many input arguments.')
% assert(nargout<3,'Too many output arguments.')

assert(isequal(size(x),size(z)),'X and Z must of same size.')
if isequal(size(SIG),[1 2]) || isequal(size(SIG),[1 3])
    nl = abs(SIG(1)); nc = abs(SIG(2));
else
    [nl,nc,~] = size(SIG);
end

% check if we have I/Q signals
isIQ = ~isreal(SIG);

%-- Check input parameters
if ischar(varargin{end})
    method = varargin{end};
    NArg = nargin-1;
else
    method = 'linear';
    NArg = nargin;
end

if NArg==4 % DASMTX(SIG,x,z,param)
    if isstruct(varargin{1})
        param = varargin{1};
        param = IgnoreCaseInFieldNames(param);
    else
        error('The structure PARAM is required.')
    end
    assert(isfield(param,'TXdelay'),...
        'A TX delay vector (PARAM.TXdelay or DELAYS) is required.')
    delaysTX = param.TXdelay;
else % NArg=5: DASMTX(SIG,x,z,delaysTX,param)
    delaysTX = varargin{1};
    param = varargin{2};
    assert(isstruct(param),'The structure PARAM is required.')
    param = IgnoreCaseInFieldNames(param);
end
if isfield(param,'TXdelay') % DASMTX(SIG,x,z,delaysTX,param)
    assert(isequal(delaysTX(:),param.TXdelay(:)),...
        'If both specified, PARAM.TXdelay and DELAYS must be equal.')
end
assert(isequal(size(x),size(z)),'X and Z must be of same size.')


%-- Interpolation method
if ~ismember(lower(method),{'nearest','linear','quadratic','lanczos3','5points','lanczos5'})
    error('METHOD must be ''nearest'', ''linear'', ''quadratic'', ''Lanczos3'', ''5points'' or ''Lanczos5''.')
end

%-- Propagation velocity (in m/s)
if ~isfield(param,'c')
    param.c = 1540; % longitudinal velocity in m/s
end

%-- Sampling frequency (in Hz)
if ~isfield(param,'fs')
    error('A sampling frequency (PARAM.fs) is required.')
end

%-- f-number
if ~isfield(param,'fnumber')
    param.fnumber = 0; % f-number (default = full aperture)
elseif isempty(param.fnumber)
    % do nothing:
    % The f-number will be determined automatically
else
    assert(isscalar(param.fnumber) && isnumeric(param.fnumber),...
        'If not empty, PARAM.fnumber must be a scalar.')
    assert(param.fnumber>=0,...
        'PARAM.fnumber must be non-negative.')
end

%-- Acquisition start time (in s)
if ~isfield(param,'t0')
    param.t0 = 0; % acquisition start time in s
end

%-- Pitch & width or kerf (in m)
if ~isfield(param,'pitch')
    error('A pitch value (PARAM.pitch) is required.')
end
if isempty(param.fnumber)
    % NOTE:
    % An element width or a kerf width is required if the f-number is []
    if isfield(param,'width') && isfield(param,'kerf')
        assert(abs(param.pitch-param.width-param.kerf)<eps,...
            'The pitch must be equal to (kerf width + element width).')
    elseif isfield(param,'kerf')
        param.width = param.pitch-param.kerf;
    elseif isfield(param,'width')
        param.kerf = param.pitch-param.width;
    else
        error(['An element width (PARAM.width) or kerf width ',...
            '(PARAM.kerf) is required if PARAM.fnumber = [].'])
    end
    ElementWidth = param.width;
end

%-- -6dB bandwidth (in %)
if ~isfield(param,'bandwidth')
    param.bandwidth = 75;
end
if isempty(param.fnumber)
    assert(param.bandwidth>0 && param.bandwidth<200,...
        'The fractional bandwidth at -6 dB (PARAM.bandwidth, in %) must be in ]0,200[')
end

%-- Radius of curvature (in m)
% for a convex array
if ~isfield(param,'radius')
    param.radius = Inf; % default = linear array
end
RadiusOfCurvature = param.radius;
isLINEAR = isinf(RadiusOfCurvature);

%-- Reception angle (in rad) -- [Advanced option for vector Doppler] --
if ~isfield(param,'RXangle') || param.RXangle==0
    isRXangle = false;
    param.RXangle = 0;
else
    %- this option is not available for convex arrays
    assert(isinf(RadiusOfCurvature),...
        'PARAM.RXangle must be 0 with a convex array.')
    %-
    isRXangle = true;
    cosRX = cos(param.RXangle);
    tanRX = tan(param.RXangle);
    if isscalar(param.RXangle)
        cosRX = cosRX*ones(size(x));
        tanRX = tanRX*ones(size(x));
    end
end

%-- Passive imaging
if ~isfield(param,'passive')
    param.passive = false;
else
    assert(islogical(param.passive),...
        'PARAM.passive must be a boolean (false or true)')
end

%-- Number of elements
assert(numel(delaysTX)==nc,...
    'DELAYS and/or PARAM.TXdelay must be vectors of length size(SIG,2).')
% Note: param.Nelements can be required in other functions of the
%       Matlab Ultrasound Toolbox
if isfield(param,'Nelements')
    assert(param.Nelements==numel(delaysTX),...
        'PARAM.TXdelay or DELAYS must be of length PARAM.Nelements.')
end

%-- Locations of the transducer elements (!if PARAM.elements is given!)
if isfield(param,'elements')
    assert(ismatrix(param.elements) &&...
        any(ismember([1 2],size(param.elements))),...
    ['PARAM.elements must be a vector, a 2-row or a 2-column matrix ',...
    'that contains the x- (and optionally z-)locations of the transducer elements.'])
    if size(param.elements,1)==2
        % param.elements is a two-row matrix
        xe = param.elements(1,:);
        ze = param.elements(2,:);
    elseif size(param.elements,2)==2
        % param.elements is a two-column matrix
        xe = param.elements(:,1)';
        ze = param.elements(:,2)';
    else
        % param.elements is a vector
        xe = param.elements(:)';
        ze = 0;
    end
    isPARAMelements = true;
else
    isPARAMelements = false;
end

%-- Center frequency (in Hz)
if isIQ
    if isfield(param,'fc')
        if isfield(param,'f0')
            assert(abs(param.fc-param.f0)<eps,...
                ['A conflict exists for the center frequency:',13,...
                'PARAM.fc and PARAM.f0 are different!'])
        end
    elseif isfield(param,'f0')
        param.fc = param.f0; % Note: param.f0 can also be used
    else
        error('A center frequency (PARAM.fc) is required with I/Q data.')
    end
    wc = 2*pi*param.fc;
end


%-------------------------------%
% end of CHECK THE INPUT SYNTAX %
%-------------------------------%



%-- Centers of the tranducer elements (x- and z-coordinates)
if isLINEAR
    % Linear array
    if ~isPARAMelements
        xe = ((0:nc-1)-(nc-1)/2)*param.pitch;
        ze = 0;
    end
else
    % Convex array
    chord = 2*RadiusOfCurvature*...
        sin(asin(param.pitch/2/RadiusOfCurvature)*(nc-1));
    h = sqrt(RadiusOfCurvature^2-chord^2/4);
    THe = linspace(atan2(-chord/2,h),atan2(chord/2,h),nc);
    if ~isPARAMelements
        ze = RadiusOfCurvature*cos(THe);
        xe = RadiusOfCurvature*sin(THe);
        % ze = ze-RadiusOfCurvature;
        ze = ze-h;
    end
end
% note: THe = angle of the normal to element #e with respect to the z-axis


% some parameters
fs = param.fs; % sampling frequency
c = param.c; % propagation velocity


% Interpolations of the TX delays for a better estimation of dTX
idx = ~isnan(delaysTX);
assert(sum(abs(diff(idx)))<3,...
    'Several simultaneous sub-apertures are not allowed.')
if ~param.passive
    nTX = nnz(idx); % number of transmitting elements
    if nTX>1
        idxi = linspace(1,nTX,4*nTX);
        xTi = interp1(xe(idx),idxi,'spline');
        if isLINEAR
            zTi = 0;
        else
            zTi = interp1(ze(idx),idxi,'spline');
        end
        delaysTXi = interp1(delaysTX(idx),idxi,'spline');
    else
        xTi = xe(idx);
        if isLINEAR, zTi = 0; else, zTi = ze(idx); end
        delaysTXi = delaysTX(idx);
    end
end


%-- f-number (determined automatically if not given)
% The f-number is determined from the element directivity
% See the paper "So you think you can DAS?"
if isempty(param.fnumber)
    lambdaMIN = c/(param.fc*(1+param.bandwidth/200));
    RXa = abs(param.RXangle);
    % Note: in Matlab, sinc(x) = sin(pi*x)/(pi*x)
    f = @(th,width,lambda)...
        abs(cos(th+RXa)*sinc(width/lambda*sin(th+RXa))-0.71);
    opt = optimset('TolX',pi/100);
    alpha = fminbnd(@(th) f(th,ElementWidth,lambdaMIN),0,pi/2-RXa,opt);
    param.fnumber = 1/2/tan(alpha);
end
fNum = param.fnumber;


t0 = param.t0;
x = x(:); z = z(:);

%----
% Migration - diffraction summation (Delay & Sum, DAS)
%----

%-- TX distances
% For a given location [x(k),z(k)], one has:
% dTX(k) = min(delaysTXi*c + sqrt((xTi-x(k)).^2 + (zTi-z(k)).^2));
% In a compact matrix form, this yields:
if ~param.passive
    dTX = min(delaysTXi*c + sqrt((xTi-x).^2 + (zTi-z).^2),[],2);
else
    dTX = 0;
end

%-- RX distances
% For a given location [x(k),z(k)], one has:
% dRX(k) = sqrt((xT-x(k)).^2 + (zT-z(k)).^2);
% In a compact matrix form, this yields:
dxT = x-xe; dzT = z-ze;
dRX = sqrt(dxT.^2 + dzT.^2);

%-- Travel times
tau = (dTX+dRX)/c;

%-- Corresponding fast-time indices
idxt = (tau-t0(:)')*fs + 1;

%-- In-range indices:
switch lower(method)
    case 'nearest', I = idxt>=1 & idxt<=nl;
    case 'linear', I = idxt>=1 & idxt<=nl-1;
    case 'quadratic', I = idxt>=1 & idxt<=nl-2;
    case 'lanczos3', I = idxt>=2 & idxt<=nl-2;
    case '5points', I = idxt>=3 & idxt<=nl-2;
    case 'lanczos5', I = idxt>=3 & idxt<=nl-3;
end

%-- Aperture (using the f-number):
if fNum>0
    if ~isRXangle
        if isfinite(RadiusOfCurvature)
            % -- for a convex array
            Iaperture = abs(asin(dxT./dRX)-THe)<=atan(1/2/fNum);
        else
            % -- for a linear array
            % For a given location [x(k),z(k)], one has:
            % Iaperture = abs(xT-x(k)) <= z(k)/2/fNum;
            % In a compact matrix form, this yields:
            Iaperture = abs(dxT)<=(z/2/fNum);
        end

    else % [Advanced option for vector Doppler: Reception angle]
        % -- ONLY for a rectilinear array
        % For a given location [x(k),z(k)], one has:
        % Iaperture = abs(xT-x(k)-z(k)*tanRX(k)) <= z(k)/cosRX(k)/2/fNum;
        % In a compact matrix form, this yields:
        Iaperture = abs(dxT-z.*tanRX(:))<=(z./cosRX(:)/2/fNum);
    end
    I = I&Iaperture;
end

% subscripts to linear indices (instead of using SUB2IND)
idx = idxt + (0:nc-1)*nl;
idx = idx(I);


%-- Let's fill in the sparse DAS matrix
switch lower(method) % Interpolation Method
    
    case 'nearest' %-- Nearest neighbor interpolation (1-point method)
        
        %-- DAS matrix
        [i,~] = find(I);
        j = round(idx);
        if isIQ % phase rotation (if I/Q signals)
            s = exp(1i*wc*tau(I));
        else
            s = 1;
        end
        if isfield(param,'TransposeDASMatrix')
            % -- DASMTX has been called by the function DAS --
            %
            % The smallest DAS matrix (in terms of memory) is returned.
            % (Matlab stores sparse matrices in compressed sparse column format)
            if numel(x)>nl*nc
                param.TransposeDASMatrix = false;
                M = sparse(i,j,s,numel(x),nl*nc);
                % M is a [numel(x)]-by-[nl*nc] sparse matrix
            else
                param.TransposeDASMatrix = true;
                M = sparse(j,i,s,nl*nc,numel(x));
                % M is a [nl*nc]-by-[numel(x)] sparse matrix
            end
        else
            M = sparse(i,j,s,numel(x),nl*nc);
            % M is a [numel(x)]-by-[nl*nc] sparse matrix
        end
        
        
    case 'linear' %-- Linear interpolation (2-point method)
        
        idxf = floor(idx);
        idx = idxf-idx;
        
        %-- DAS matrix
        [i,~] = find(I);
        j = [idxf; idxf+1];
        s = [idx+1; -idx];
        if isIQ % phase rotation (if I/Q signals)
            tau = tau(I);
            s = s.*exp(1i*wc*[tau;tau]);
        end
        if isfield(param,'TransposeDASMatrix')
            % -- DASMTX has been called by the function DAS --
            %
            % The smallest DAS matrix (in terms of memory) is returned.
            % (Matlab stores sparse matrices in compressed sparse column format)
            if numel(x)>nl*nc
                param.TransposeDASMatrix = false;
                M = sparse([i;i],j,s,numel(x),nl*nc);
                % M is a [numel(x)]-by-[nl*nc] sparse matrix
            else
                param.TransposeDASMatrix = true;
                M = sparse(j,[i;i],s,nl*nc,numel(x));
                % M is a [nl*nc]-by-[numel(x)] sparse matrix
            end
        else
            M = sparse([i;i],j,s,numel(x),nl*nc);
            % M is a [numel(x)]-by-[nl*nc] sparse matrix
        end
        
        
    case 'quadratic' %-- quadratic interpolation (3-point method)
        
        idxf = floor(idx);
        idx = idx-idxf;
        
        %-- DAS matrix
        [i,~] = find(I);
        j = [idxf; idxf+1; idxf+2];
        s = [(idx-1).*(idx-2)/2;
            -idx.*(idx-2);
            idx.*(idx-1)/2];
        if isIQ % phase rotation (if I/Q signals)
            tau = tau(I);
            s = s.*exp(1i*wc*repmat(tau,3,1));
        end
        if isfield(param,'TransposeDASMatrix')
            % -- DASMTX has been called by the function DAS --
            %
            % The smallest DAS matrix (in terms of memory) is returned.
            % (Matlab stores sparse matrices in compressed sparse column format)
            if numel(x)>nl*nc
                param.TransposeDASMatrix = false;
                M = sparse(repmat(i,3,1),j,s,numel(x),nl*nc);
                % M is a [numel(x)]-by-[nl*nc] sparse matrix
            else
                param.TransposeDASMatrix = true;
                M = sparse(j,repmat(i,3,1),s,nl*nc,numel(x));
                % M is a [nl*nc]-by-[numel(x)] sparse matrix
            end
        else
            M = sparse(repmat(i,3,1),j,s,numel(x),nl*nc);
            % M is a [numel(x)]-by-[nl*nc] sparse matrix
        end
        
        
    case 'lanczos3' %-- 3-lobe Lanczos interpolation (4-point method)
        
        idxf = floor(idx);
        idx = idx-idxf;
        
        %-- DAS matrix
        [i,~] = find(I);
        j = [idxf-1; idxf; idxf+1; idxf+2];
        s = [sinc(idx+1).*sinc((idx+1)/2);
            sinc(idx).*sinc(idx/2);
            sinc(idx-1).*sinc((idx-1)/2);
            sinc(idx-2).*sinc((idx-2)/2)];
        if isIQ % phase rotation (if I/Q signals)
            tau = tau(I);
            s = s.*exp(1i*wc*repmat(tau,4,1));
        end
        if isfield(param,'TransposeDASMatrix')
            % -- DASMTX has been called by the function DAS --
            %
            % The smallest DAS matrix (in terms of memory) is returned.
            % (Matlab stores sparse matrices in compressed sparse column format)
            if numel(x)>nl*nc
                param.TransposeDASMatrix = false;
                M = sparse(repmat(i,4,1),j,s,numel(x),nl*nc);
                % M is a [numel(x)]-by-[nl*nc] sparse matrix
            else
                param.TransposeDASMatrix = true;
                M = sparse(j,repmat(i,4,1),s,nl*nc,numel(x));
                % M is a [nl*nc]-by-[numel(x)] sparse matrix
            end
        else
            M = sparse(repmat(i,4,1),j,s,numel(x),nl*nc);
            % M is a [numel(x)]-by-[nl*nc] sparse matrix
        end
        
        
    case '5points' %-- 5-point least-squares parabolic interpolation
        
        idxf = floor(idx);
        idx = idx-idxf;
        
        %-- DAS matrix
        [i,~] = find(I);
        j = [idxf-2; idxf-1; idxf; idxf+1; idxf+2];
        idx2 = idx.^2;
        s = [1/7*idx2-1/5*idx-3/35;
            -1/14*idx2-1/10*idx+12/35;
            -1/7*idx2+17/35;
            -1/14*idx2+1/10*idx+12/35;
            1/7*idx2+1/5*idx-3/35];
        if isIQ % phase rotation (if I/Q signals)
            tau = tau(I);
            s = s.*exp(1i*wc*repmat(tau,5,1));
        end
        if isfield(param,'TransposeDASMatrix')
            % -- DASMTX has been called by the function DAS --
            %
            % The smallest DAS matrix (in terms of memory) is returned.
            % (Matlab stores sparse matrices in compressed sparse column format)
            if numel(x)>nl*nc
                param.TransposeDASMatrix = false;
                M = sparse(repmat(i,5,1),j,s,numel(x),nl*nc);
                % M is a [numel(x)]-by-[nl*nc] sparse matrix
            else
                param.TransposeDASMatrix = true;
                M = sparse(j,repmat(i,5,1),s,nl*nc,numel(x));
                % M is a [nl*nc]-by-[numel(x)] sparse matrix
            end
        else
            M = sparse(repmat(i,5,1),j,s,numel(x),nl*nc);
            % M is a [numel(x)]-by-[nl*nc] sparse matrix
        end
        
        
    case 'lanczos5' %-- 5-lobe Lanczos interpolation (6-point method)
        
        idxf = floor(idx);
        idx = idx-idxf;
        
        %-- DAS matrix
        [i,~] = find(I);
        j = [idxf-2; idxf-1; idxf; idxf+1; idxf+2; idxf+3];
        s = [sinc(idx+2).*sinc((idx+2)/2);
            sinc(idx+1).*sinc((idx+1)/2);
            sinc(idx).*sinc(idx/2);
            sinc(idx-1).*sinc((idx-1)/2);
            sinc(idx-2).*sinc((idx-2)/2);
            sinc(idx-3).*sinc((idx-3)/2)];
        if isIQ % phase rotation (if I/Q signals)
            tau = tau(I);
            s = s.*exp(1i*wc*repmat(tau,6,1));
        end
        if isfield(param,'TransposeDASMatrix')
            % -- DASMTX has been called by the function DAS --
            %
            % The smallest DAS matrix (in terms of memory) is returned.
            % (Matlab stores sparse matrices in compressed sparse column format)
            if numel(x)>nl*nc
                param.TransposeDASMatrix = false;
                M = sparse(repmat(i,6,1),j,s,numel(x),nl*nc);
                % M is a [numel(x)]-by-[nl*nc] sparse matrix
            else
                param.TransposeDASMatrix = true;
                M = sparse(j,repmat(i,6,1),s,nl*nc,numel(x));
                % M is a [nl*nc]-by-[numel(x)] sparse matrix
            end
        else
            M = sparse(repmat(i,6,1),j,s,numel(x),nl*nc);
            % M is a [numel(x)]-by-[nl*nc] sparse matrix
        end

        
end

end


function structArray = IgnoreCaseInFieldNames(structArray)

switch inputname(1)
    case 'param'
        fieldLIST = {'attenuation','baffle','bandwidth','c','fc',...
            'fnumber','focus','fs','height','kerf','movie','Nelements',...
            'passive','pitch','radius','RXangle','RXdelay',...
            'TXapodization','TXfreqsweep','TXnow','t0','width'};
    case 'options'
        if isstruct(structArray)
            fieldLIST = {'dBThresh','ElementSplitting',...
                'FullFrequencyDirectivity','FrequencyStep','ParPool',...
                'WaitBar'};
        else
            return
        end
end

OldFieldNames = fieldnames(structArray);
tmp = lower(OldFieldNames);
assert(length(tmp)==length(unique(tmp)),...
    ['The structure ' upper(inputname(1)),...
    ' contains duplicate field names (when ignoring case).'])

[idx,loc] = ismember(lower(fieldLIST),tmp);
idx = find(idx); loc = loc(idx);
for k = 1:length(idx)
    tmp = eval(['structArray.' OldFieldNames{loc(k)}]); %#ok
    structArray = rmfield(structArray,OldFieldNames{loc(k)});
    eval(['structArray.' fieldLIST{idx(k)} ' = tmp;'])
end

end

function [M,param] = RunTheExample

%-- Generate RF signals using a phased-array transducer

% Phased-array @ 2.7 MHz:
param = getparam('P4-2v');

% TX time delays (80-degree-wide diverging wave)
dels = txdelay(param,0,80/180*pi);

% Scatterers' position:
xs = [(-1:0.5:1)*4e-2 zeros(1,5)];
zs = [ones(1,5)*6e-2 (2:2:10)*1e-2];

% Backscattering coefficient
BSC = [ones(1,9) 0];

% RF signals:
param.fs = 4*param.fc; % sampling frequency
RF = simus(xs,zs,BSC,dels,param);

% Plot the RF signals
figure
subplot(121)
plot((RF(:,1:7:64)/max(RF(:))+(1:10)*2)',...
    (0:size(RF,1)-1)/param.fs*1e6,'k')
set(gca,'XTick',(1:10)*2,'XTickLabel',int2str((1:7:64)'))
title('RF signals')
xlabel('Element number'), ylabel('time (\mus)')
xlim([0 22]), axis ij


%-- Demodulation and beamforming

% Demodulation
IQ = rf2iq(RF,param);

% Beamforming grid
[th,r] = meshgrid(linspace(-40,40,128)/180*pi+pi/2,...
    linspace(1,9,256)*1e-2);
[x,z] = pol2cart(th,r);

% DAS matrix and beamformed IQ
M = dasmtx(1i*size(IQ),x,z,dels,param);
IQb = reshape(M*IQ(:),size(x));

% Beamformed image
subplot(122)
pcolor(x*100,z*100,abs(IQb).^.5)
colormap(gray)
shading interp, axis equal ij tight
title('Gamma-compressed image')
xlabel('[cm]'), ylabel('[cm]')
hold on
% position of the scatterers (in cm)
plot(xs*100,zs*100,'ro')

end

